% 飞行器参数
% 创建与 X 相同长度的时间向量
time = linspace(0, simulation_time, length(X));
mass = 100; % 千克
drag_coefficient = 0.5;
lift_coefficient = 0.5;
wing_area = 10; % 平方米
gravity = 9.81; % 米/秒^2
density = 1.225; % 千克/立方米（标准海平面的空气密度）

% 初始条件
initial_x = 0; % 米
initial_y = 0; % 米
initial_vx = 100; % 米/秒
initial_vy = 0; % 米/秒
initial_theta = 0; % 弧度

% 仿真参数
time_step = 0.1; % 秒
simulation_time = 100; % 秒

% 初始化状态变量
x = initial_x;
y = initial_y;
vx = initial_vx;
vy = initial_vy;
theta = initial_theta;

% 初始化数据存储数组
X = zeros(1, simulation_time/time_step + 1);
Y = zeros(1, simulation_time/time_step + 1);
VX = zeros(1, simulation_time/time_step + 1);
VY = zeros(1, simulation_time/time_step + 1);
THETA = zeros(1, simulation_time/time_step + 1);

% 仿真循环
for i = 1:(simulation_time/time_step + 1)
    t = (i - 1) * time_step;
    
    % 计算力
    lift_force = 0.5 * lift_coefficient * density * vx^2 * wing_area;
    drag_force = 0.5 * drag_coefficient * density * vx^2 * wing_area;
    weight_force = mass * gravity;
    
    % 计算加速度
    ax = (lift_force * sin(theta) - drag_force * cos(theta)) / mass;
    ay = -(lift_force * cos(theta) + drag_force * sin(theta) + weight_force) / mass;
    
    % 更新状态变量
    x = x + vx * time_step;
    y = y + vy * time_step;
    vx = vx + ax * time_step;
    vy = vy + ay * time_step;
    theta = theta + (vy / vx) * time_step;
    
    % 记录数据
    X(i) = x;
    Y(i) = y;
    VX(i) = vx;
    VY(i) = vy;
    THETA(i) = theta;
end

% 绘制曲线图
figure;
subplot(3, 1, 1);
plot(time, X, 'b-');
xlabel('Time (seconds)');
ylabel('Horizontal Position (meters)');
title('Horizontal Position vs. Time');

subplot(3, 1, 2);
plot(time, Y, 'b-');
xlabel('Time (seconds)');
ylabel('Vertical Position (meters)');
title('Vertical Position vs. Time');

subplot(3, 1, 3);
plot(time, VX, 'b-', time, VY, 'r-');
xlabel('Time (seconds)');
ylabel('Velocity (meters/second)');
title('Velocity vs. Time');
legend('Horizontal Velocity', 'Vertical Velocity');